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PREFACE 


This  Memorandum  is  a  result  of  RAND's  continuing  study  of 
guidance  and  orbit  mechanics.  It  was  written  to  elaborate  on  the 
analytical  details  of  the  computational  equations  for  a  Kalman  filter, 
discussed  in  RM-A241-PR,  Tlie  Application  of  Dynamic  Programming  to 
Satellite  Intercept  and  Rendezvous  Problems.  F.  T.  Smith,  August  1964. 
Such  a  filter  was  used  to  process  observational  data  and  to  estimate 
relative  state  variables  for  the  optimal  control  process. 

In  addition,  this  Memorandum  discusses  the  application  of  the 
Kalman  filter  to  orbit  determination  and  orbit  correction  processes 
when  these  processes  arc  described  in  terms  of  two-body  orbital 
parameters.  It  may  be  of  use  in  studies  of  the  application  of 
dynamic  programming  to  flight  control  problems. 
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SUMMARY 

This  Memorandum  applies  dynamic  programming  to  the  derivation  of 
the  computational  equations  for  a  Kalman  filter. 

Ihe  problem  of  using  a  set  of  ill>conditioned  observations  is 
considered.  For  example,  when  observations  of  a  satellite  taken  from 
a  single  ground  station  are  spaced  too  close  together  in  time,  the 
accuracy  of  the  determination  of  the  satellite's  orbital  parameters 
is  poor.  Such  a  set  of  observations  is  said  to  be  ill-conditioned. 

The  general  set  of  computational  equations  derived  is  interpreted 
in  terms  of  an  orbit  determination  problem,  an  orbit  correction  pro¬ 
cess,  and  a  satellite  rendezvous  problem. 
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I.  INTRODUCTION 

This  Memorandum  uses  the  techniques  of  dynamic  programming  to 
formulate  the  Kalman  filtering  problem  and  to  obtain  the  necessary 
computational  equations.  The  discussion  is  less  general  than  that 
in  Ref.  1,  and  the  point  of  view  here  is  slanted  more  toward  astro* 
dynamical  applications.  The  application  of  dynamic  programming  to 
this  type  of  problem  is  discussed  in  Refs.  2  and  3.  This  Memorandum 
goes  into  much  more  detail  concerning  applications  than  do  Refs.  2 
and  3. 

The  method  of  least  squares  is  briefly  discussed  to  establish 

notation.  The  problem  is  then  formulated  as  a  multi-stage  decision 

process,  and  dynamic  progrontnlng  concepts  are  applied  to  obtain  the 

necessary  computational  equations.  A  method  of  treating  ill-condi- 
(4) 

tinned  systems  is  incorporated  into  the  computational  equations. 
The  application  of  the  estimation  process  to  the  determination  of  a 
set  of  two-body  orbital  parameters  from  a  sequence  of  observations  is 
then  discussed,  l.e.,  the  orbit  determination  problem.  The  discus¬ 
sion  is  then  extended  to  feedback  control  problems  where  it  is  de¬ 
sired  to  estimate  time-varying  corrections  to  two-body  orbital  pa¬ 
rameters  and  system  state  vectors. 
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II.  THE  METHOD  OF  LEAST  SQUARES 
Consider  a  set  of  m  observations  of  one  or  more  physical  quan¬ 
tities  which  have  been  corrupted  by  noise.  If  y^  Is  the  true  value 
of  each  quantity  at  the  time  of  the  1-th  observation,  and  v^  Is  the 
error  due  to  noise,  then  the  observed  value  Is 

Zi  «  Vi  +  yi,  1  -  1,  • •  •  ,  m 


Assume  that  the  y^  are  related  to  a  set  of  n  parameters  x 
ues  are  to  be  determined  from  the  relation 


J 


whose  val- 


* 

where 


y  »  Hx 


m  X  1  vector  )  y^^, 
n  X  1  vector  fx^^. 


H  B  m  X  n  matrix  whose  elements  may  be 
time  varying 

Since  the  true  values  of  the  observed  quantities  are  corrupted  by 
noise,  the  best  we  can  hope  to  do  Is  to  solve  the  equation 


z  ■  Hx  +  V 


This  equation  may  be  written  as 


where 


V  «  z  -  Hx 


The  superscript  T  denotes  the  transpose  of  a  vector  or  matrix. 


3 


We  shall  assume  that  vector  v  represents  a  gausslan  random  pro> 
cess  with  zero  mean.  According  to  the  principle  of  least  squares  the 
best  estimate  of  vector  x  will  be  that  for  which  the  quadratic  form 

T  -1 

J  -  V  R  V 

Is  a  minimum,  where  R  Is  the  covariance  matrix  of  the  gausslan  process. 
The  well-known  solution  to  this  problem  Is 

X  -  hV^z 

If  this  solution  Is  substituted  Into  the  quadratic  fora  representing 
J  we  obtain 

Min  J  -  fz  -  H(hV^H)‘^  hV^z]  r'^[z  -  H(hV^)  hV^z] 

=  H(H^r"^H)“^  H^r’^]  z  -  z’^Q(m)z 

If  matrix  H  Is  square  and  non-singular,  Min  J  Is  zero. 

Suppose  that  we  now  acquire  a  new  obsexrvation  z  and  we  wish 

nrt*! 

to  use  z  ,  to  improve  further  our  estimate  of  x.  We  can  form  an  mfl 
mf  1  ^ 

system  of  equations  and  solve  as  above  to  obtain  x  where  H  Is  now  an 
(nri-1)  X  n  matrix. 

As  we  combine  more  and  more  observations  the  estimate  x  approaches 
closer  to  the  true  vector  x,  but  the  computational  process  becomes 
more  and  more  Inefficient  since  the  dimensions  of  H  and  R  grow  with  each 
observation.  It  is  desirable,  therefore,  to  find  n  more  efficient  meth¬ 
od  of  utilizing  new  observations. 
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III.  THE  APPLICATION  OF  DYNAMIC  PROGRAMMING 


The  process  of  estimating  the  components  of  vector  x  becomes  a 
multi-stage  decision  process  when  new  observations  are  combined  se¬ 
quentially  to  obtain  successively  improved  estimates  of  x.  Further, 

each  new  decision  requires  a  new  minimization  operation. 

•k 

We  shall  call  the  vector 


z 


the  observation  vector.  At  each  stage  of  the  process  vector  z  in¬ 
creases  its  dimension  by  one  due  to  the  addition  of  a  new  observation. 
The  earlier  components  of  z,  of  course,  remain  unchanged.  At  each 
stage  of  the  process  the  decision  involves  choosing  x  such  that 


where 


J  =  v^(m)  R  ^(m)  v(m) 
m 

T 

VW  -  [vj,  ....  vj 


is  minimized  where  m  is  the  number  of  the  observations  involved. 

It  is  clear  that  the  estimate  of  x  depends  on  the  number  and  the 
values  of  the  components  of  z(m)  ,  as  will  the  value  of  J^.  Instead 


To  clarify  notation  regarding  the  dimension  of  vectors  z  and  v, 
we  shall  let  z(m)  and  v(m)  denote  the  fact  that  z  and  v  are  m-dimen- 
sional  vectors.  Vector  x  is  assumed  to  be  always  of  dimension  n. 

•kic 

The  notation  R(m)  implies  an  m  x  m  matrix. 
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of  referring  to  an  N-stage  decision  process  we  shall  refer  to  an  m- 

observation  estimation  process  and  define  the  minimum  value  of  J  as 

m 


The  cost  of  an  m-observat ion  estimation  process 
based  on  an  m>dimensional  observation  vector 
z(m)  with  an  optimal  estimation  policy  being 
used. 


By  an  optimum  estimation  policy  is  meant  the  choice  of  successive 
estimates  of  x  such  that  the  performance  index 


J  -  v^m) 


R  ^(m)  v(m) 


m  ■  n ,  n+1 ,  • • • 


is  minimized.  Hie  reason  for  index  m  starting  at  m  ■  n  is  that  a 
minimum  of  n  observations  are  required  to  obtain  a  useful  estimate 


of  the  n  components  of  vector  x.  It  is  assumed  that  the  components 
of  v(m)  are  independent,  i.e.,  R  ^(m)  is  a  m  x  m  diagonal  matrix. 


Let  a  new  observation  be  added  as  the  (m+l)-th  component  to 

vector  z(m).  A  new  component  ^be  uncertainty  associated  with 

observation  z^^^,  is  also  added  to  vector  v(m).  If  R^^  ^  is  the  var- 

'fc 

lance  associated  with  v  then  we  can  write  J  as 

m+1  m+1 

T  -1 

*^m+l  '^m+1  ^m+1  '^m+l  "^m 

-  '’l+l  Vrt  ''■n+l  +  ’‘’‘("J 


*  -1 
Since  v  and  R  are 
m+1  m+1 

The  vector-matrix  notation  is 


scalars  equals 

retained  for  convenience  since  it  is  pos¬ 


sible  to  add  several  observations  at  a  time  to  improve  the  estimate, 

and  R  would  then  be  a  matrix. 
m+1 
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The  estimate  of  x  that  minimizes  J  denoted  by  (x)  will  be 

m+i  m-n+1 

* 

different  from  that  minimizing  J  .  Let  the  new  estimate  minimizing 

fit 

J  . ,  be  expressed  as  (x)  +  (^x)  .  V/hen  the  components  of  vector 

m+1  m-n  '  'm-n 

(x)  are  changed  to  (x)  +  Ax  a  new  value  for  J  >  f  [z(m)]  is 

m-n  m-n  mm 

obtained 

T 

^m  '  -  H(m)  +  AxJ}  R‘^(m)  {zCm) 

■  “("J  +  '^]] 


Exp2tnding  this  quadratic  form  and  properly  grouping  terms  yields 


J  *  2*^(0!)  R  ^(m)  z(m)  -  2  z^(m)  R  ^(m)  H(m)  (x) 
m  m-n 

n  HC™)  (x)_  _ 

m-n  m-n 


+  2  Ax^  R”^(m)  H(m) 

-  H^(m)  R  ^(m)  z^(m)J  +  Ax^  H^(m)  R  ^(m)  H(m)  Ax 


The  first  three  terms  clearly  are  f 

m 

vanishes  from  the  definition  of  (x) 

m-n 


z(m)  J  and  the  bracketed  term 
(See  Section  II).  This  leaves 


J 

m 


T  T 

Ax  H  (m) 


R"^(m)  H(m)  Ax  +  f 


m 


The  notation  (x)  denotes  the  (m-n)-th  estimate  of  vector  x 
m-n 

and  is  used  to  distinguish  the  (m-n)-th  estimate  from  the  (m-n)-th 
component  denoted  by  x^ 
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and  J  becomes 
m+1 


m+1 


T 

{*«+!  -  Vl  +  ^*]}  Cl  {'m+1  '  «m+l  +  '«]} 


+  Ax^  H^(ra)  R  ^(m)  H(m)  Ax  +  f ^  |^z(ni)  J 

We  can  now  apply  the  principle  of  optimality  to  derive  a  recurrence 
relation  for  f^  |^z(m)J.  For  this  problem  the  principle  of  optimality 
can  be  stated  as  follows: 

An  optimal  sequence  of  estimations  of  the  parameter 
vector,  (x)j^,  (x)2»  has  the  property  that  what¬ 

ever  choice  is  made  for  (x)j^,  the  remaining  sequence 
of  estimates  (x)2,  (x)^*  *•*»  must  constitute  an  op¬ 
timal  sequence  with  regard  to  the  observations  from 
which  (x)^  is  estimated. 

We  obtain  the  following  recurrence  relations: 

*m+l  [^<"'+1)]  '  "te  {'m+1  -  "m+l  [<'>m.n  +  '^]}  Cl  {'m+1 

”  ^m+1  H(m)  Ax 


+  f 


m 


Tne  symbol  H 


vation  z  ... 

m+1 

is  a  matrix. 


If  z 


ra+1 

m+1 


denotes  the  n  x  1  vector  associated  with  obser- 

denotes  several  simultaneous  observations  H 

m+1 
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where 


Q(n)  z(n) 


m  ■  n ,  n+1 ,  • • • 

The  first  term  inside  the  braces  represents  the  direct  cost  of  using 
the  (m-f-l)-th  observation.  The  second  cost  represents  the  change  in 
the  cost  for  the  previous  m  observations  due  to  changing  the  estimate 
of  X. 

Noting  that  f^  j^z(m)J  is  independent  of  Ax,  we  minimize  the  terms 

inside  the  braces  by  differentiating  with  respect  to  the  components  of 

Ax.  Setting  the  result  equal  to  zero  and  solving  for  (Ax)  yields 

m-n 


«<")  *  Vl]’*  "Ih  [ 


'm+1 


-  H  U.1  (x)  1 

ra+1  '  'm-nj 


The  new  estimate  for  x  is  then 


(x)  , ,  *  (x)  +  (Ax) 

'  m-n+1  m-n  m-n 


It  is  clear  that  can  be  expressed  as 


m+1 


1 

o 

B 

1 

i 

v(m) 

1 

+ 

O 

J 

''m+l 

v^(m+l)  R  ^(m+1)  v(m+l) 
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where  R  ^m+l)  denotes  the  partitioned  matrix. 
Since 


v(m) 

z(m) 

H(m) 

'^m+l 

^m+1 

»m+l 

it  is  clear  from  Section  II  that  f^^j^[z(m+l)  ]  must  be  given  by 


f^_^l|^2(m+l)  J  «  z^(m+l)  Q(m+1)  z(m+l) 


where 


.-1. 


Q(nifl)  -  R  (m+1) 

-  R  ^(m+1)  H(m+1)  |^H'^(ra+l)  R  ^(m+1)  H(m+l)J  H^(m+1)  R  ^(m+1) 


H(m+1) 


H(ra) 


H 


m+1 


It  follows  from  the  nature  of  as  expressed  by  the  equation 
J_  *  H^(ra)  R  ^(m)  H(m)  Ax  +  f _  |^z(m)J 


m 


m 


that  each  correction  to  the  estimate  of  vector  x  for  each  stage  of  the 
decision  process  can  be  computed  by  the  equation 


«<'">  ^  'tl  ■‘mil  "m+l]  '  Cl  [ 


m+1 


-  H  (x)  1 

m+1  m-nj 


10 


where 


(x)  . ,  «  (x)  +  (Ax) 

'  ra-n+1  '  m-n  'm-n 


m  ■  n ,  n-fl ,  •  • 


To  avoid  a  matrix  inversion  process  we  follow  the  notation  of 
Ref.  1  and  define 


P"^  »  H^(m)  R‘^(m)  H(ra) 

m 


-1  -1  T  -1 

P  ^  -  P  ^  H 

m+1  m  m+1  m+1  m+1 


Then  by  the  lemma  given  in  Ref,  1  we  have 


P  =  P  -  P  P  +  R  H  P 

m+1  m  m  m+1  m+1  m  m+1  m+1  m+1  m 


If  we  now  replace  our  index  notation  by  k  »  m-n  we  have  the  following 
set  of  recurrence  relations: 


<^>k+l  •  <^>k 

^^^k  ‘  ^k+1  \+l  \+l  [^k+1  "  “k+1 

\+l  ^  \  \  ”k+l  ^\+l  \  \+l  \+l^  \+l  ^k 


k  «  0,  1,  2, 


This  can  be  proved  by  direct  substitution  of  the  expressions 

Cl  ‘‘m+1  Cl  ‘‘m+l  - 
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In  order  to  use  these  recurrence  relations  it  is  necessary  to  have 
a  set  of  initial  conditions  (x)^  and  as  well  as  R^,  •••. 

The  quantity  can  represent  a  single  observation  or  a  vec¬ 

tor  whose  components  represent  several  observations.  The  computa¬ 
tional  advantage  of  utilizing  observations  one  at  a  time  is  that 
T 

^\+l  \  \+l  \+l^  \+l  are  scalars  and  matrix  inversion  is 

eliminated  entirely. 


The  selection  of  suitable  initial  conditions  will  be  discussed  in 


the  sections  to  follow. 
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IV.  ILL -CONDITIONED  SYSTPIS 

Consider  the  equation  relating  the  true  observation  vector  y  and 
the  parameter  vector  x  given  in  Section  II. 

y  »  H  X 

In  order  that  this  equation  have  a  unique  solution  the  rank  of  matrix 
H  must  equal  the  dimension  of  vector  x.  Equivalently,  we  may  require 
that  the  determinant  of  matrix  i.e.,  the  Gram  determinant  of  H,  be 

non-zero.  While  this  condition  satisfies  the  theoretical  requirement 
for  a  unique  solution,  it  may  not  satisfy  computational  requirements 
for  finding  such  a  solution. 

Vftien  the  Gram  determinant  has  a  sufficiently  small  numerical  value, 
the  computation  of  the  components  of  vector  x  is  subject  to  considerable 

9 

error  and  the  set  of  observations  defining  vector  y  is  said  to  be  ill- 
conditioned. 

The  recurrence  equations  given  in  the  last  section  may  converge 
too  rlowly,  or  not  at  all,  if  they  are  applied  to  an  ill-conditioned 
set  of  observations. 

A  method  for  treating  ill-conditioned  linear  systems  is  given  by 

(it) 

Bellman,  Kalaba,  and  Lockett.  The  recurrence  relations  given  in 
Section  III  can  be  modified  to  incorporate  one  of  these  computational 
techniques  for  ill-conditioned  systems. 

Consider  first  the  one-stage  process  in  Section  II.  We  modify 
the  performance  index  to  be 

-|T  ^ 

■  l^z(m)  -  H(ra)xJ  R  (m)  |^z(m)  -  H(m)xJ 

+  X(x-c)^  (x-c) 

where  c  is  an  n  x  1  vector  chosen  so  that  the  norm  | |x-c  ]  ]  is  small. 
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If  we  expand  J  ,  differentiate  with  respect  to  the  components  of  vec- 
m 

tor  X,  set  the  result  equal  to  zero,  and  solve  for  x,  we  obtain 


A 

X 


^-1 


j^H^(m)  R“\m)  H(m)  +  XlJ  H^(m)  R‘^(m)  z(m) 


+  X  R  ^(m)  H(m)  +  u] 


as  the  best  least  squares  estimate.  Reference  4  describes  a  method  of 
successive  approximations  where  vector  c  is  taken  as  (x)^^  and  x  be¬ 
comes  (x)2.  We  then  have  the  recurrence  equation 


^*^k+l  “  H(m)  +  Xl]  {H^(ra)  R’^(m)  z(m)  +  >.(x)j^]' 


which,  as  shown  in  Ref.  4,  converges  to  x. 

Assume  that  a  new  observation  z  becomes  available  whose  un- 

m+1 

certainty  is  We  know  that  the  optimal  estimate  of  x  will  change. 

The  value  of  J  for  (x)  +  Ax  becomes 

m  m-n 


T 

Jjj,  -  {z(n>)  -  H(m)  +  AxJ|  R"^(m)  |z(m)  -  H(m)  +  -]} 


+  [<*>n,-n  +  4*  -  ']  [(*)„.„  +  to  -  c] 


If  the  subscripts  and  arguments  are  temporarily  dropped  for 
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notatlonal  convenience  we  can  expand  J  as 

m 


T  — 

=  j^z  -  H(x  +  R”^rz  -  H(x  +  Ax)J  +  X(x  +  Ax  -  c)^  (x  +  Ax  -  c) 


zV^z  -  2(x  +  Ax)  hV^z  +  (x  +  Ax)^  (hV^H  +  XI)  (x  +  Ax) 


-  2X(x  +  Ax)^  c  +  Xc^c 


We  next  expand  the  terms  involving  x  +  Ax  and  regroup  the  terms 
to  take  advantage  of 

[hV)  R"^(m)  H(m)  +  Xl]  -  H^(m)  R‘\m)  z(m)  -  Xc  -  0 


We  obtain 


J  -  zV^z  -  2x^V^z  +  x^(H^R‘^H  +  XI)  X  -  2AxVr'^s 


+  2Ax'*^(H^r”^H  +  XI)  X  +  Ax^(H^R"^H  +  XI)  Ax  -  2XAx^c 


„.aT  ,  ,  T 
2Xx  c  +  Xc  c 


T_-l  aT„T_-1  ^  aT  r/wTo-l„  ^  *  itTo-l  ^  1 

zR  z-xHR  z+x  I  (HR  H  +  XI)  x  -  H  R  z  -  Xc 


+  2Ax^  +  XI)  X  -  H^r'^z  -  Xc] 


+  Ax^(H^r"^H  +  Xl)  Ax  -  Xx^c  +  Xc^c 
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Hie  third  and  fourth  vanish.  We  next  substitute  for  x  in 

the  second  term  and  obtain 


T  -1  T  -1  T  T  -1  -1  T  -1 

=  z  R  z  -  (H  R  z  +  Xc)  (H  R  H  +  XI)  H  R  z 


+  Ax^(H^R‘^H  +  XI)  Ax  -  Xx^c  +  Xc'^c 


Expanding  the  second  term  gives 


J  =  z 
m 


^  R'^  -  R"^K(hV^H  +  XI)  h’^R'^J  z  -  Xc^(hV^H  +  XI)  hV^z 


+  Ax^(hV^H  +  XI)  Ax  -  Xx^c  +  Xc^c 


2  T  T  *1  *1  Ta 

We  next  add  and  subtract  the  term  X  c  (H  R  H  +  Xl)  c  -  Xc  x  and 
obtain  for  J 


J  *  z 
m 


j^R“^  -  R"^H(hV^H  +  XI) H^R’^] 


Xc^  [(hV^H  +  XI) H^R“^z  +  X(hV4  +  XI) c  -  x] 


-  X^c^(hV^H  +  XI)  c  +  Xc’x  -  Xx^c  +  Xc^c 


+  Ax^(hV^H  +  XI)  Ax 
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Hie  second  term  vanishes  and  two  terms  cancel  leaving 

J  =  Tr"^  -  R“^H(h\‘^H  +  XI)  hV^I  z 

m  ^  J 

+  Xc^  [^I  -  X(H^R"^H  +  XI) c 
+  Ax^(H^r"^H  +  XI)  Ax 
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Comparing  this  recurrence  relation  with  that  in  Section  III 

for  |^z(m+l)J  shows  that  the  terns  involving  vector  Ax  in  the 

T  -1 

two  expressions  differ  only  by  having  H  (m)  R  (m)  H(m)  replaced 
T  -1 

by  H  (m)  R  (m)  H(m)  +  Xl.  If  we  make  a  similar  substitution  in 

the  equation  for  (Ax)  we  obtain 

m-n 


«<"•>  +  “^1  Vl  +  “m+l  [= 


-  H  (x)  1 

m+1  '  'm-nj 


Define 


H^(in)  R"^(m)  H(m)  +  XI 


-1  T  -1 

p  +  H* R  H  . , 
m  m'fl  m+1  m+1 


From  the  lemma  in  Ref.  1,  we  have,  as  in  Section  III 


P  -  P  -  P  (H  P  +  R  H  P 

m+1  m  m  m+1  m+l  m  m+1  m+l  m+1  m 


We  can  thus  use  the  final  set  of  recurrence  relations  given  in 
Section  III  for  ill-conditioned  systems.  The  difference  for  ill- 
conditioned  systems  concerns  the  initial  matrix  P^,  which  is  given 


P^  »  H^^Cn)  R’^(n)  H(n)  +  XI 
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and  the  Initial  estimate  (x)  »  which  is  obtained  by  the  method  of 

o 

successive  approximations  as  described  in  Ref.  4  using  the  first  n 
observations. 
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V.  THE  ORBIT  DETERMINATION  PROBLEM 
■Hie  determination  of  the  definitive  orbit  of  a  celestial  object 
usually  involves  a  differential  correction  process  which  relates  ob¬ 
servation  residuals  to  corrections  of  a  set  of  two-body  orbital  pa- 
•k 

rameters.  Tliis  problem  can  be  quite  easily  treated  as  a  Kalman  fil¬ 
tering  problem,  since  the  observational  data  are  usually  large  in 
number  and  the  sources  can  be  separated  in  space  and  time.  This  ap¬ 
plication  was  originally  proposed  by  Swerling.^^^ 

As  an  example,  consider  the  differential  correction  process 
discussed  in  Ref.  7  .  The  residuals  of  a  sequence  of  azimuth  and 
elevation  angle  observations  are  to  be  related  to  the  corrections  to 
I  >  .  of  six  two-body  orbital  parameters.  These  corrections  can  be 
represented  as  an  orbital  parameter  correction  vector 


^  T 

Ap  =  [Aa/a,  A(e  cos  E^) ,  A(e  sin  E^) ,  Au^^,  Av^,  Aw^] 

u  (8) 

where 


a  is  the  semi-major  axis 


e  is  the  eccentricity 


E  is  the  value  of  the  eccentric  anomaly  at  t  ■  t 
o  o 

Au  ,  Av  ,  Aw  are  small  rotations  about  a  set  of  unit  vectors 
o  o’  o 

i?  ,  V  ,  W  which  define  the  orientation  of  the 
o  o 

orbit  in  space. 

The  small  rotations  Au  ,  Av  ,  and  Aw  can  be  related  to  small  changes 


in  the  unit  vectors  AU  ,  A^  ,  AW.  The  vector  Ap  corresponds  to  vector 

o  o 

X  in  Sections  II  and  III. 


See  Chapter  11  of  Ref.  b  for  a  discussion  of  this  problem. 
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The  true  values  of  the  observation  residuals  are  represented  as 
the  ra  X  1  vector 


Aq  -  [Aq^,  Aq2,  Aq^] 


where 


Aq.  *  -  p.  AA.  cos  h. 
ill  i 


i  =  1,  3,  *•*,  m-l 


'“'i+i  *  Pi  ^*>1 


The  quantities  and  Ah^  are  the  true  values  of  the  observation 
residuals  for  the  azimuth  and  elevation  angles. 

Within  the  limits  of  linearity  the  vectors  Aq  and  Ap  are  re¬ 
lated  by  the  relation 


Aq  ■  CAp 


where  C  is  an  m  x  6  matrix  whose  elements  are  functions  of  time  and 
the  elements  of  the  preliminary  two-body  orbit  being  corrected.  The 
explicit  expressions  giving  the  elements  of  matrix  C  arc  rather  in¬ 
volved,  Since  they  arc  derived  in  detail  in  Ref.  7,  they  will  not  be 
given  hen. 

If  we  assume  that  the  observations  are  corrupted  by  white  gaussian 
noise  of  zero  mean,  then  the  above  vector  matrix  equation  becomes 

V  ■  z  -  CAp 

where  the  observation  vector  in  Sections  II  and  III  is  given  by 

z  «  V  +  Aq 


V 


[V 


1’ 


^m^ 

m 


T 


and 
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The  ccxnponents  of  obsexvation  vector  z  are 


\  -  *cl>  '■i 


1^  1|  3|  ID*  1 


z, *  p. (h  -  h  ,) 
1+1  1  oi  Cl 


where  A  . ,  h  ,  are  the  observed  and  A  . ,  h  ,  are  the  computed  values 

Ol  01  Cl  Cl 

of  azimuth  and  elevation  angles  at  times  t..  Angles  A  .  and  h  .  are 

1  Cl  Cl 

based  on  the  values  of  the  orbital  parameters  of  the  preliminary  or¬ 
bit  being  corrected,  as  are  the  factors  and  cos  h^. 

By  associating  APj^  with  x^,  6pj^  with  and  matrix  C  with  ma¬ 

trix  H  the  recurrence  relations  in  Section  III  become 


=  ‘■k+i  ^+1  \li  (Vi  - 


^Pk+i  •  *Pk  ^  ‘Pk 


k  =  1 ,  2,  ... 


^k+1  “  ^k+1  ^^k+1  \  ^k+1  \+l^  ^k+1 


where  the  subscript  k  refers  to  the  k-th  estimate. 

If  we  assume  that  each  observation  results  in  a  pair  of  residuals 
in  azimuth  and  elevation,  then  matrix  is  2  x  6 ,  matrix  is  6  x  6, 
and  matrix  is  2  x  2.  Vector  will  be  2  x  1. 

If  we  take  the  first  three  pairs  of  observations  we  can  form  the 
equation  r  -i  r  -i 


Ap  *  CAp 
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where 


T 

[-  ■  *cl)  -  h^pj  .  I  .  1,2,3 


We  can  compute  the  initial  estimate  of  Ap  as 


-  C  ^  z 


Azimuth  and  elevation  angle  observations  are  independent.  Assume 

2 

they  have  equal  variances  o  which  are  constant  for  all  observations. 
Then 


R,  -  R2  -  R3 


0 


I, 


Let 


R  = 


0 

0 


0 

0 


I. 


Then  since  is  the  covariance  matrix  for  it  can  be  computed  from 

2  T  -1 
=  a  (C^  C) 

It  is  also  possible  to  let  P^^  be  a  diagonal  matrix  whose  elements 
are  estimated  on  the  basis  of  experience,  and  to  let  Ap^^  be  the  null 
vector. 


If  the  matrix  C  is  singular  then  the  observations  are  not  in¬ 
dependent  and  the  whole  process  breaks  down.  The  process  discussed 
in  Section  IV  can  be  used  if  the  observations  arc  ill-conditioned. 
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VI.  THE  FEEDBACK  CONTROL  PROBLEM 

THE  ORBIT  TRANSFER  PROBLEM  AS  AN  OPTIMAL  REGULATOR  PROBLEM 

An  orbit  transfer  process  can  be  represented  as  a  feedback  con¬ 
trol  process  and  treated  as  an  optimal  rei,ulator  problem.  Consider 
that  the  space  vehicle  is  in  an  orbit  defined  by  p^  at  time  t^  and 

it  is  desired  to  transfer  it  to  an  orbit  defiled  by  p  at  time  t  +  t. 

T  o 

The  symbols  p^  and  p^  represent  vectors  having  as  components  the  sets 

of  values  of  the  two-body  orbital  parameters  describing  the  osculating 

orbits  to  the  vehicle's  trajectories  at  t  and  t  +  t.  We  define  an 

o  o 

orbital  cirrection  vector  at  time  t  as 

Ap(ty  =  -  P(t) 

where  p(t)  is  the  vector  representing  the  values  of  the  osculating 
orbital  parameters  at  time  t.  The  object  of  the  orbit  transfer  pro¬ 
cess  is  then  to  drive  the  orbital  parameter  correction  vector  to  zero 

at  time  t  +  T.  By  some  optimization  process  the  optimal  control  vec- 
o 

tor  can  be  expressed  as  a  functi  jn  of  the  orbital  parameter  correction 

vector  Ap(t).  This  closes  a  feedback  loop  and  expresses  the  orbit 

transfer  problem  as  an  optimal  regulator  problem.  It  also  implies  the 

necessity  for  computing  the  components  of  Ap(t)  as  functions  of  certain 

observed  quantities,  since  the  components  of  Ap(t)  are  not  directly 
ic 

observable. 

The  orbit  transfer  problem  can  be  formulated  in  phase  space.  In 

this  case  vectors  p^  and  p^  are  replaced  by  s  and  s^.  The  components 

of  vectors  s  and  s^  are  the  rectangular  components  of  position  and 

velocity  in  the  initial  and  terminal  trajectories.  A  relative  system 
* 

Assuming  Ap(t)  is  inaccurately  known. 
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state  vector  can  then  be  defined  as 
As(t)  «  s^(t)  -  s(t) 

*  [Ax(t),  Ax(t),  Ay(t),  Ay(t),  Az(t),  A2(t)]^ 

In  phase  space  the  values  of  the  components  of  the  system  state  vectors 

s  and  s^  change  with  time,  with  or  without  thrust  being  applied  to  the 

vehicle.  The  object  of  the  orbit  transfer  problem  is  then  to  make  state 

vector  s  coincide  with  vector  s  at  time  t  +  t. 

r  o 

T5ie  vector  s^(t)  can  represent  the  motion  of  a  hypothetical  point 
mass  moving  along  the  desired  terminal  trajectory  of  the  orbit  transfer 
process,  or  it  can  represent  the  motion  of  an  actual  space  vehicle  that 
is  to  be  intercepted  or  rendezvoused  with.  In  the  intercept  case  it 
may  be  desired  only  to  drive  the  relative  position  components  Ax(t), 
Ay(t),  and  A2(t)  to  zero  at  time  t^  +  t. 

As  the  components  of  As(t)  are  usually  not  directly  observable, 
it  is  necessary  to  compute  them  from  observable  quantities  in  order 
that  the  optimal  control  process  can  be  carried  out. 

Actually,  in  phase  space,  the  orbit  transfer  is  more  like  a  servo 
problem  than  a  regulator  problem  since  As(t)  refers  to  a  moving  origin. 
However,  under  certain  circumstances  it  can  be  treated  approximately 
as  a  regulator  problem  by  neglecting  certain  small  terms  or  defining 
As(t)  as  As(t)  -  t)  -  s(t). 

In  the  remainder  of  this  section  the  application  of  Kalman  fil¬ 
tering  theory  to  the  estimation  of  orbital  parameter  correction  vec¬ 
tors  Ap(t)  and  relative  system  state  vectors  As(t)  will  be  discussed. 
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THE  ESTIMATION  OF  ORBITAL  PARAMETER  CORRECTION  VECTORS 

In  the  discussion  to  follow  it  is  assumed  that  the  durations  of 
time  involved  are  short  enough  so  that  the  space  vehicle  essentially 
follows  a  two-body  orbit  while  no  thrust  is  being  applied,  i.e.,  the 
effect  of  perturbations  from  natural  causes  is  assumed  negligible. 

Under  these  conditons  the  motion  of  the  body  is  defined  by  a  set  of 
six  constant  orbital  parameters.  The  rectangular  position  and  velocity 
of  the  space  vehicle  can  be  determined  at  any  particular  time  from  these 
orbital  parameters  by  applying  suitable  transformation  equations.  How¬ 
ever,  when  thrust  is  applied  to  the  space  vehicle  the  values  of  these 
orbital  parameters  change  continuously  with  time.  At  each  instant  of 
time  the  resulting  set  of  orbital  parameters  defines  an  osculating  conic. 
This  conic  represents  the  path  in  space  along  which  the  vehicle  would 
move  if  the  thrust  were  to  vanish.  An  orbit  transfer  can  thus  be  con¬ 
sidered  as  the  rime  behavior  of  a  two-body  orbit  that  is  continuously 
deformed  in  size,  shape,  and  orientation  in  space  from  some  initial  con- 
fiouration  to  some  terminal  configuration.  Because  of  the  nonlinearity 
of  the  set  of  differential  equatic  s  defining  the  time  behavior  of  the 
orbital  parameters,  the  estimation  problem  will  be  treated  as  a  discrete 
process  and  linear  approximations  used. 

The  orbit  determination  problem  involves  obtaining  the  estimated 
values  of  a  set  of  six  constant  parameters  from  a  sequence  of  observa¬ 
tions.  The  feedback  control  problem  involves  determining  the  estimated 
values  of  a  set  of  six  time-varying  orbital  parameters  at  discrete  in¬ 
stants  of  time  during  the  orbit  transfer  process.  Thus,  some  relation 
is  necessary  for  determining  how  the  orbital  parameters  change  their 
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values  with  time. 

A  piece-wise  linear  solution  to  an  optimal  control  problem  formu¬ 
lated  in  terms  of  orbital  parameters  and  obtained  by  dynamic  program¬ 
ming  is  given  in  Ref.  9.  The  linear  transformation  relating  the  or¬ 
bital  parameter  correction  vector  at  time  t^^  to  that  at  time  is 

given  by 

^^k+l  “  \  •••,  N-1 

where  is  a  6  x  3  matrix  whose  elements  are  evaluated  from  the  com¬ 
ponents  of  the  osculating  orbital  parameter  vector  at  tj^.  The  vector 
is  the  incremental  change  in  the  space  vehicle's  velocity  vector 
from  tj^  to  due  to  the  application  of  thrust.  From  Ref.  9 

where  a  6  x  6  matrix,  is  determined  from 

\  ”  \+l  ■  \+l\  \+l  \  \  \+l 

Si'  ^ 

Substituting  from  in  the  transformation  equation  yields 

''Pk+i  \  \+i  \  *  \+i]  '"’k 


Multiplying  both  sides  by  matrix  and  solving  for  Apj^  gives 


.-1 


''■’k  ■  \  Vi  '^Pk+i 


This  relation  defines  how  the  orbit  parameter  correction  vector  changes 
from  one  stage  of  the  process  to  the  next  during  an  optimal  orbit  trans¬ 
fer. 

It  will  be  assumed  that  at  each  time  tj^  an  observation  is 

obtained,  an  orbital  parameter  correction  vector  estimate  Ap(tj^)  will 
be  computed  (Fig.  1).  The  optimal  control  process  will  utilize  this 
estimate  to  compute  a  change  in  the  vehicle's  velocity  vector.  The 
orbital  parameter  correction  vector  and  the  true  values  of  the  obser¬ 
vations  at  t,  are  related  by 
k 

Aq(tj^)  =  C(p^;  t^)  Ap(tj^) 

Matrix  C(p  ;  t,  )  is  the  m  x  6  matrix  discussed  in  Section  IV  except 

T  fC 

that  its  elements  are  evaluated  in  terms  of  the  components  of  p  .  Ihis 

T 

means  that  the  linear  approximations  involved  in  deriving  the  elements 
of  matrix  C  become  more  exact  as  the  termination  of  the  process  is  ap¬ 
proached. 

Associating  Aq(tj^^j^)  with  ^^m’  ^P^*Tc+l^ 

X  ,  and  C(p  :  t,  with  H  we  can  rewrite  the  recurrence  rela- 
m+1  k+1  m+1’ 

tions  with  a  necessary  slight  modification  as 

®Pk+l  ■  ^k+1  ®k+l  \+l  ■  ‘^k+1  \+l  \ 

VP  ■  \  ^p<'k>  ‘P<VP 

Vi  -  \  \  \  Vl  ^  VP''  Vl  \ 


where 


Rela 
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The  slight  modification  mentioned  above  involves  the  extrapolation  of 
Ap(t^)  ahead  to  time  by  using  the  Q  matrices.  The  sequence  of  Q 

matrices  is  generated  in  the  process  of  solving  the  optimal  control 
problem. 

THE  ESTIMATION  OF  RELATIVE  SYSTEM  STATE  VECTORS 

The  linearized  time-varying  differential  equation  for  the  relative 
system  state  vector  is  given^^^^  by 

^  *=  A(t)  As  -  [A(t)  -  A^]  s^(t)  -  Bu(t),  As(t^)  -  Ac 

where  A(t)  is  a  time-varying  matrix,  A^  is  a  constant  matrix  for  the 
target  in  a  circular  orbit,  B  is  a  constant  matrix,  and  u(t)  is  the 
control  vector.  For  a  quadratic  system  performance  index  the  optimal 
control  vector  will  be  a  linear  function  of  the  relative  system  state 
vector 

u(t)  -  ^  ^(y(T  +  t^-  t)  As(t)  +  e(T  +  t^-  t)  s^(t)J 

Substituting  this  expression  for  u(t)  in  the  differential  equation 
above  yields 

^  =  F(t  +  t  -  t)  As  +  G(t  +  t  -  t)  s  (t) 
dt  o  o  r 

F(t  +  t^-  t)  -  A(t)  -  I”  B  B^Qr(T  +  t^-  t) 

G(t  +  t  -  t)  -  A  -  A(t)  -  Y  B  B^P(T  +  t  -  t) 
o  r  X  o 


where 
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The  solution  of  this  time-varying  differential  equation  is  given 


by 

t 

As(t)  -  $(t.  t^)  As(t^)  +  j  $(t,  §)  G(§)  s^CO  d§ 

t 

o 


*  4(t,  t^)  As(t^)  +  g(t,  t^) 


where  matrix  $(t,  t^)  satisfies  the  matrix  differential  equation 


dt 


F(t  +  t^-  t)  $  , 


$(t  ,  t  )  *  I 
^  0  o 


Ihis  solution  defines  the  time  behavior  of  the  relative  state 
vector  as  the  vehicle  moves  along  a  trajectory  determined  by  the  op¬ 
timal  control  vector. 

The  relative  system  state  variables  are  expressed  with  respect 
to  a  non-rotating  coordinate  system  and  are  not  directly  observable. 
T!ieir  estimated  values  must  therefore  be  based  on  physical  quantities 
that  can  be  observed.  There  is.  obviously  more  than  one  way  of  handling 
this  problem.  Hie  method  to  be  discussed  consists  of  two  separate  com¬ 
putational  processes. 

°  Observations  of  relative  range,  relative  range  rate, 
two  angles,  and  two  angular  rates  associated  with  the 
line  of  sight  from  s(t)  to  Sj.(t)  are  made  at  discrete 
instants  of  time.  T!ic  corresponding  values  of  the 
State  variables  are  computed  by  suitable  transforma¬ 
tion  equations. 

o  Hie  sequence  of  sets  of  computed  values  of  the  state 
variables  is  operated  on  by  a  Kalman  filter  to  op¬ 
timally  estimate  values  of  the  state  variables  up¬ 
dated  in  time  from  the  observation  times. 
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Assume  that  As(tj^),  k  »  1,  2,  •••,  ra  represents  the  computed  rela¬ 
tive  system  state  vector  based  on  observations  at  t^^,  k  »  1,  2,  •••,  m. 

The  problem  is  to  estimate  the  As(T)  where  T  >  t^^,  k  *  1,  2,  •••,  ra. 

The  relative  system  state  vector  at  t  «  T  is  related  to  the  rela¬ 
tive  system  state  vector  at  t^^  by  the  solution  of  the  differential 

equation 

As(T)  -  $(T,  tj^)  As(tj^)  + 

«  $(T,  t^)  As(t^)  +  g(T,  t^) 
or 

As(tj^)  -  $'\t.  tj^)  [as(T)  -  g(T,  tj^)] 

By  making  the  proper  associations  of  vectors  and  matrices  we  can  re¬ 
write  the  recurrence  relations  in  Section  III  as 

[6s(I)]j^  ■  ’’k+1  *  ^\+l’ 

^+1  ■  '■k  -  '■k  r»<‘k+i-  ^ 


J 

/ 


$(T.  1)  G(§)  s^(§)  d§ 
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where 

T 

g(T.  t^)  .  j  $(T.  §)  G(§)  s^(§) 
k*0  ,  1,  •••jin-l 

There  will  be  some  uncertainty  associated  with  g(T,  tj^)  due  to 
uncertainty  in  knowledge  of  the  target's  orbit.  However,  the  matrix 
G(§)  will  have  elements  that  are  either  relatively  small  or  zero,  and 
since  g(T,  tj^)  is  in  the  nature  of  a  correction  the  effects  of  un¬ 
certainties  associated  with  g(T,  tj^)  will  be  ignored. 

The  initial  quantities  required  by  the  recurrence  relations  are 
[As(T)]^  ,  ,  and  R^(tj^).  The  initial  estimate  of  [As(T)]^  can  be 

computed  from 

[As(T)]^  -  $(T,  t^)  As(t^) 


The  matrix  can  be  taken  as  a  diagonal  matrix  whose  non-zero  elements 
are  the  variances  associated  with  the  initial  values  of  the  relative 


system  state  variables. 

If  the  matrix  R  (t,  )  is  the  covariance  matrix  associated  with  the 

q  k 

set  of  observations  taken  at  time  t^  and  the  matrix  Jacobian 

matrix  associated  with  the  equations  transforming  from  the  observed 


quantities  to  the  computed  relative  system  state  variables,  then  when 
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the  observation  uncertainties  are  sufficiently  small 

R^(t^)  .  M(t^)  R^(t^)  h\) 

where  R^(t^)  is  the  covariance  matrix  for  the  relative  system  state 

variables  at  t,  . 

k 

The  process  just  described  estimates  the  values  of  the  relative 
system  state  variables  at  time  T  >  t^,  k  ■  1,  2,  •••,  m  based  on  m  sets 
of  observations.  The  feedback  control  process  requires  that  the  rela¬ 
tive  system  state  variables  be  estimated  either  continuously  or  at  a 
discrete  number  of  times  throughout  the  orbit  transfer  process.  Since 
the  computational  processes  involved  will  usually  be  performed  by  a 
digital  computer,  both  cases  are  actually  discrete,  the  difference 
being  due  to  the  frequency  with  which  As(t)  is  estimated. 

The  period  of  time  during  which  observational  data  arc  processed 
to  obtain  the  estimate  As(T),  i.e.,  0  S  t  i  T  will  be  referred  to  as 
an  estimation  interval.  It  is  assumed  that  the  orbit  trai’sfer  process 
consists  of  N  estimation  intervals,  i.e.,  As(t)  is  estimated  at  times 
T,  2T,  •••,  (N-l)T  with  m  sets  of  observations  being  made  during  each 
estimation  interval.  The  process  involved  in  estimating  As(T)  is  then 
repeated  during  each  estimation  interval,  and  it  remains  to  consider 
how  the  various  estimation  intervals  are  linked  together  so  that  all 
preceding  observational  data  arc  used  for  each  successive  estimate. 

The  initial  estimate  of  As(2T)  can  be  computed  from 


[As(2T)]^  =  K2T,  T)  [As(T)]^_^  +  g(2T,  T) 
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where 

2T 

g(2T,  T)  -  j  J(2I,  5)  G(|)  s^(5)  d5 
T 

The  initial  P  matrix  for  the  second  estimation  Interval  is  com¬ 
puted  from 

P  -  P  ,  -  P  ,  [P  ,  +  R  P  - 
m  ra-1  m-1  m-1  m**  m-1 

since 

Ht,  T)  -  $(T,  T)  -  I 
m 

The  computation  of  R  (t)  depends  only  on  the  time  of  the  associ''.ted 

s 

set  of  observations,  and  the  equation  given  above  is  used  for  all  es- 

★ 

timation  intervals. 

Imbedding  the  one-stage  estimation  process  in  the  N-stage  esti¬ 
mation  process  provides  estimates  of  the  relative  system  state  vari¬ 
ables  at  (N-1)  Instants  of  time  during  the  orbit  transfer.  When  a 
continuous  control  process  is  involved,  the  digital  computation  pro¬ 
cess  used  to  compute  As(t)  is  based  on  a  simple  linear  extrapolation 
during  each  estimation  interval;  i.e., 

As(t)  -  $(t,  nT)  As(nT)  +  g(t,  nT) ,  n  -  1 ,  2,  •  •  • ,  N-1 
where  ^ 

g(t,  nT)  .  j  »(t,  ?)  G(5)  s^(5)  dC 
nT 

★ 

Since  the  matrix  M(t)  becomes  singular  as  s(t)  approaches  s^(t) 

the  computation  process  for  R  (t)  may  have  o  be  modified  during  the 

s 

last  few  estimation  ii  ’’ervals. 
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